import numpy
cimport numpy

cimport cython

cdef extern from "math.h":
	double exp(double)
	double sqrt(double)

@cython.boundscheck(False)
def logistic0(numpy.ndarray[numpy.float64_t,ndim=1] vec, double s):
	cdef int N = vec.shape[0]

	cdef numpy.ndarray[numpy.float64_t,ndim=1] vec1
	vec1 = vec/s

	cdef numpy.ndarray[numpy.float64_t,ndim=1] probs
	probs = numpy.exp(-vec1)
	probs /= 1.0 + numpy.exp(-vec1)
	probs /= 1.0 + numpy.exp(-vec1)
	probs /= s

	return probs

def logistic0_scalar(double val,s):
	cdef double val1 = val/s
	cdef double prob = numpy.exp(-val1)
	prob /= 1.0 + exp(-val1)
	prob /= 1.0 + exp(-val1)
	prob /= s

	return prob


@cython.boundscheck(False)
def logistic_weights(numpy.ndarray[numpy.float64_t,ndim=1] vec, double coeff):
	cdef int N = vec.shape[0]
#	cdef double stdev = vec.std()
	cdef double stdev = numpy.sqrt((vec*vec).sum()/N)

	cdef numpy.ndarray[numpy.float64_t,ndim=1] u1
	u1 = numpy.zeros((N,),numpy.float)
	u1 = logistic0(vec,stdev)

	cdef numpy.ndarray[numpy.float64_t,ndim=1] u2
	u2 = numpy.zeros((N,),numpy.float)
	u2 = logistic0(vec,coeff*stdev)

	cdef numpy.ndarray[numpy.float64_t,ndim=1] weights
	weights = u1/(u1 + u2)

	cdef double p1 = logistic0_scalar(0.0,stdev)
	cdef double p2 = logistic0_scalar(0.0,coeff*stdev)

	weights *= (p1 + p2)/p1

	return weights



@cython.boundscheck(False)
def redescend_coefficients(numpy.ndarray[numpy.float64_t,ndim=1] vec,
			double coeff):

	cdef numpy.ndarray[numpy.float64_t,ndim=1] weights
	weights = logistic_weights(vec,coeff)

	return vec*weights


#@cython.boundscheck(False)
#def outlier_contrast(numpy.ndarray[numpy.float64_t,ndim=2] obs):
#	cdef Py_ssize_t N = obs.shape[0]
#	cdef Py_ssize_t P = obs.shape[1]
#
#	cdef double w0,we,wc
#
#	cdef int i,j
#
#	cdef numpy.ndarray[numpy.float64_t,ndim=2] contrast
#	contrast = numpy.zeros((obs.shape[0],obs.shape[1]),numpy.float)
#
#	cdef numpy.ndarray[numpy.float64_t,ndim=2] middle
#	middle = numpy.zeros((3,3),numpy.float)
#	middle[0,0] = middle[0,2] = middle[2,0] = middle[2,2] = -1.0/numpy.sqrt(2.0)
#	middle[0,1] = middle[2,1] = middle[1,0] = middle[1,2] = -1.0
#
#	middle[1,1] = -middle.sum()
#	middle /= numpy.sqrt((middle*middle).sum())
#
#	cdef numpy.ndarray[numpy.float64_t,ndim=2] corner
#	corner = numpy.zeros((2,2),numpy.float)
#	corner[1,1] = -1.0/numpy.sqrt(2.0)
#	corner[0,1] = corner[1,0] = -1.0
#
#	corner[0,0] = -corner.sum()
#	corner /= numpy.sqrt((corner*corner).sum())
#
#	cdef numpy.ndarray[numpy.float64_t,ndim=2] edge
#	edge = numpy.zeros((2,3),numpy.float)
#	edge[1,0] = edge[1,2] = -1.0/numpy.sqrt(2.0)
#	edge[0,0] = edge[0,2] = edge[1,1] = -1.0
#
#	edge[0,1] = -edge.sum()
#	edge /= numpy.sqrt((edge*edge).sum())
#
#
#	w0 = corner[0,0]
#	we = corner[0,1]
#	wc = corner[1,1]
#
#	# Top left
#	contrast[0,0] = w0*obs[0,0] + we*(obs[1,0] + obs[0,1]) + wc*obs[1,1]
#
#	# Top right
#	contrast[0,P-1] = w0*obs[0,P-1] + wc*obs[1,P-2]
#	contrast[0,P-1] += we*(obs[1,P-1] + obs[0,P-2])
#
#	# Bottom left
#	contrast[N-1,0] = w0*obs[N-1,0] + wc*obs[N-2,1]
#	contrast[N-1,0] += we*(obs[N-2,0] + obs[N-1,1])
#
#	# Bottom right
#	contrast[N-1,P-1] = w0*obs[N-1,P-1] + wc*obs[N-2,P-2]
#	contrast[N-1,P-1] += we*(obs[N-1,P-2] + obs[N-2,P-1])
#
#	w0 = edge[0,1]
#	we = edge[0,0]
#	wc = edge[1,0]
#
#	# Top and bottom
#	for j in range(1,P-1):
#		contrast[0,j] = w0*obs[0,j] + wc*(obs[1,j-1] + obs[1,j+1])
#		contrast[0,j] += we*(obs[0,j-1] + obs[0,j+1] + obs[1,j])
#
#		contrast[N-1,j] = w0*obs[N-1,j] + wc*(obs[N-2,j-1] + obs[N-2,j+1])
#		contrast[N-1,j] += we*(obs[N-1,j-1] + obs[N-1,j+1] + obs[N-2,j])
#
#
#	for i in range(1,N-1):
#		contrast[i,0] = w0*obs[i,0] + wc*(obs[i-1,1] + obs[i+1,1])
#		contrast[i,0] += we*(obs[i-1,0] + obs[i+1,0] + obs[i,1])
#
#		contrast[i,P-1] = w0*obs[i,P-1] + wc*(obs[i-1,P-2] + obs[i+1,P-2])
#		contrast[i,P-1] += we*(obs[i-1,P-1] + obs[i+1,P-1] + obs[i,P-2])
#
#	w0 = middle[1,1]
#	we = middle[0,1]
#	wc = middle[0,0]
#
#	for i in range(1,N-1):
#		for j in range(1,P-1):
#			contrast[i,j] = w0*obs[i,j]
#			contrast[i,j] += we*(obs[i,j-1] + obs[i,j+1])
#			contrast[i,j] += we*(obs[i-1,j] + obs[i+1,j])
#			contrast[i,j] += wc*(obs[i-1,j-1] + obs[i-1,j+1])
#			contrast[i,j] += wc*(obs[i+1,j-1] + obs[i+1,j+1])
#
#
#	return contrast
#
#
#@cython.boundscheck(False)
#def outlier_contrast_weighted(numpy.ndarray[numpy.float64_t,ndim=2] obs,
#			numpy.ndarray[numpy.float64_t,ndim=2] wts):
#	cdef int N = obs.shape[0]
#	cdef int P = obs.shape[1]
#
#	cdef double cwt = 1.0/numpy.sqrt(2.0)
#	cdef double val, the_sum, the_sumsq
#
#	cdef numpy.ndarray[numpy.float64_t,ndim=2] contrast
#	contrast = numpy.zeros((obs.shape[0],obs.shape[1]),numpy.float)
#
#	cdef numpy.ndarray[numpy.float64_t,ndim=2] C
#
#	C = numpy.zeros((2,2),numpy.float)
#
#	cdef int i, j
#
#
#	# Top left
#	C.fill(0.0)
#	C[0,1] = -wts[0,1]
#	C[1,0] = -wts[1,0]
#	C[1,1] = -cwt*wts[1,1]
#	C[0,0] = -C.sum()
#	C /= numpy.sqrt((C*C).sum())
#
#	contrast[0,0] = C[0,0]*obs[0,0] + C[0,1]*obs[0,1]
#	contrast[0,0] += C[1,0]*obs[1,0] + C[1,1]*obs[1,1]
#
#
#	# Top right
#	C.fill(0.0)
#	C[0,0] = -wts[0,P-2]
#	C[1,1] = -wts[1,P-1]
#	C[1,0] = -cwt*wts[1,P-2]
#	C[0,1] = -C.sum()
#	C /= numpy.sqrt((C*C).sum())
#
#	contrast[0,P-1] = C[0,0]*obs[0,P-2] + C[0,1]*obs[0,P-1]
#	contrast[0,P-1] += C[1,0]*obs[1,P-2] + C[1,1]*obs[1,P-1]
#
#	# Bottom left
#	C.fill(0.0)
#	C[0,0] = -wts[N-2,0]
#	C[1,1] = -wts[N-1,1]
#	C[0,1] = -cwt*wts[N-2,1]
#	C[1,0] = -C.sum()
#	C /= numpy.sqrt((C*C).sum())
#
#	contrast[N-1,0] = C[0,0]*obs[N-2,0] + C[0,1]*obs[N-2,1]
#	contrast[N-1,0] += C[1,0]*obs[N-1,0] + C[1,1]*obs[N-1,1]
#
#	# Bottom right
#	C.fill(0.0)
#	C[0,1] = -wts[N-2,P-1]
#	C[1,0] = -wts[N-1,P-2]
#	C[0,0] = -cwt*wts[N-2,P-2]
#	C[1,1] = -C.sum()
#	C /= numpy.sqrt((C*C).sum())
#
#	contrast[N-1,P-1] = C[0,0]*obs[N-2,P-2] + C[0,1]*obs[N-2,P-1]
#	contrast[N-1,P-1] += C[1,0]*obs[N-1,P-2] + C[1,1]*obs[N-1,P-1]
#
#
#	C = numpy.zeros((2,3),numpy.float)
#
#	# Top 
#	for j in range(1,P-1):
#		C.fill(0.0)
#		C[0,0] = -wts[0,j-1]
#		C[0,2] = -wts[0,j+1]
#		C[1,1] = -wts[1,j]
#		C[1,0] = -cwt*wts[1,j-1]
#		C[1,2] = -cwt*wts[1,j+1]
#		C[0,1] = -C.sum()
#		C /= numpy.sqrt((C*C).sum())
#
#		contrast[0,j] = (C*obs[0:2,j-1:j+2]).sum()
#
#	# Bottom
#	for j in range(1,P-1):
#		C.fill(0.0)
#		C[1,0] = -wts[N-1,j-1]
#		C[1,2] = -wts[N-1,j+1]
#		C[0,1] = -wts[N-2,j]
#		C[0,0] = -cwt*wts[N-2,j-1]
#		C[0,2] = -cwt*wts[N-2,j+1]
#		C[1,1] = -C.sum()
#		C /= numpy.sqrt((C*C).sum())
#
#		contrast[N-1,j] = (C*obs[N-2:N,j-1:j+2]).sum()
#
#
#	C = numpy.zeros((3,2),numpy.float)
#
#	# Left
#	for i in range(1,N-1):
#		C.fill(0.0)
#		C[0,0] = -wts[i-1,0]
#		C[2,0] = -wts[i+1,0]
#		C[1,1] = -wts[i,1]
#		C[0,1] = -cwt*wts[i-1,1]
#		C[2,1] = -cwt*wts[i+1,1]
#		C[1,0] = -C.sum()
#		C /= numpy.sqrt((C*C).sum())
#
#		contrast[i,0] = (C*obs[i-1:i+2,0:2]).sum()
#
#	# Right
#	for i in range(1,N-1):
#		C.fill(0.0)
#		C[0,1] = -wts[i-1,P-1]
#		C[2,1] = -wts[i+1,P-1]
#		C[1,0] = -wts[i,P-2]
#		C[0,0] = -cwt*wts[i-1,P-2]
#		C[2,0] = -cwt*wts[i+1,P-2]
#		C[1,1] = -C.sum()
#		C /= numpy.sqrt((C*C).sum())
#
#		contrast[i,P-1] = (C*obs[i-1:i+2,P-2:P]).sum()
#
#
#	C = numpy.zeros((3,3),numpy.float)
#
#
##herehere
#	# Middle
#	for i in range(1,N-1):
#		for j in range(1,P-1):
##			C.fill(0.0)
#
#			C[1,1] = 0.0
#			C[0,1] = -wts[i-1,j]
#			C[1,0] = -wts[i,j-1]
#			C[1,2] = -wts[i,j+1]
#			C[2,1] = -wts[i+1,j]
#			C[0,0] = -cwt*wts[i-1,j-1]
#			C[0,2] = -cwt*wts[i-1,j+1]
#			C[2,0] = -cwt*wts[i+1,j-1]
#			C[2,2] = -cwt*wts[i-1,j+1]
#			C[1,1] = -C.sum()
#			C /= sqrt((C*C).sum())
#
##			contrast[i,j] = (C*obs[i-1:i+2,j-1:j+2]).sum()
#			val = C[0,0]*obs[i-1,j-1] + C[0,1]*obs[i-1,j] + C[0,2]*obs[i-1,j+1]
#			val += C[1,0]*obs[i,j-1] + C[1,1]*obs[i,j] + C[1,2]*obs[i,j+1]
#			val += C[2,0]*obs[i+1,j-1] + C[2,1]*obs[i+1,j] + C[2,2]*obs[i+1,j+1]
#			contrast[i,j] = val
#
#
#	return contrast
